/*
 * Copyright 2003 Program of Computer Graphics, Cornell University
 *     580 Rhodes Hall
 *     Cornell University
 *     Ithaca NY 14853
 * Web: http://www.graphics.cornell.edu/
 * 
 * Not for commercial use. Do not redistribute without permission.
 */
#include "stdafx.h"

#include "Disk.h"

PrimitiveList *Disk::compile(int num,const GeometryObject **array,PoolAllocator *pool) const
{
//_asm int 3; //force breakpoint
  PackedDiskList *list = PackedDiskList::poolAlloc(num,pool);
  for(int k=0;k<num;k++) {
    assert(array[k]->getType()==GEOM_DISK);
    const Disk *t = (const Disk*)array[k];
    list->setDisk(k,t->centerX,t->centerY,t->centerZ,t->normalX,t->normalY,t->normalZ,
      t->uX,t->uY,t->uZ,t->radius,t->id());
  }
  return list;
}

void Disk::expandBoxToInclude(AxisAlignedBox *box) const 
{
  F32vec4 norm = getNormal();
  F32vec4 cent = getCenter();
  F32vec4 delta = F32vec4(1.0f) - norm*norm;
  delta = F32vec4(radius)*sqrt(simd_max(delta,F32vec4(0.0f)));
  box->addPoint(cent-delta);
  box->addPoint(cent+delta);
}

AxisAlignedBox Disk::getBoundingBox() const 
{
  F32vec4 norm = getNormal();
  F32vec4 cent = getCenter();
  F32vec4 delta = F32vec4(1.0f) - norm*norm;
  delta = F32vec4(radius)*sqrt(simd_max(delta,F32vec4(0.0f)));
  return AxisAlignedBox(cent-delta,cent+delta);
}

bool Disk::intersectsRay(const Ray &ray, IntersectionRecord *hitRec) const
{
  F32vec4 norm = getNormal();
  F32vec4 cent = getCenter();
  F32vec4 uAxis = getUAxis();
  //first form the denominator
  F32vec1 denom = rcp_nr(dotProduct(norm,ray.direction));
  //next form the v-axis
  F32vec4 vAxis = crossProduct(norm,uAxis);
  //S is center minus the ray origin
  F32vec4 S = cent - ray.origin;
  //t is the parametric parameter along the ray of the intersection
  //t = (S dot C) / (D dot C)
  F32vec1 t = dotProduct(S,norm)*denom;
  //B is S X D
  F32vec4 B = crossProduct(S,ray.direction);
  //u is the first barycentric coordinate of the intersection
  //u = (vAxis dot B) / (D dot C)
  F32vec1 u = dotProduct(vAxis,B)*denom;
  //v is the second barycentric coorinate of the intersection
  //v = -(uAxis dot B) / (D dot C)
  F32vec1 negV = dotProduct(uAxis,B)*denom;
  //test to see if intersection is valid
  F32vec1 rad2Test = u*u + negV*negV;
  F32vec1 radius2( radius*radius );
  F32vec1 test = cmple(rad2Test,radius2)&cmpge(t,F32vec1(ray.tStart))&cmple(t,F32vec1(ray.tEnd));
  //I haven't thought of a better way of extracting mask from F32vec1 so we convert it to
  //a F32vec4, move the mask to an integer and then ignore everything except the low bit
  if (move_mask(F32vec4(test))&0x01) {
    F32vec1 invRadius = rcp_nr(F32vec1(radius));
    u *= invRadius;
    negV *= invRadius;
    Store_F32vec1(&hitRec->t,t);
    hitRec->primitiveId = id();
    Store_F32vec1(&hitRec->u,u);
    Store_F32vec1(&hitRec->v,F32vec1(0)-negV);
    hitRec->w = 0;
    return true;
  }
  return false;
}

bool Disk::intersectsBox(const AxisAlignedBox &box) const
{
  //trivial accept if center is in the box
  if (box.containsPoint(getCenter())) return true;

  //we use a simple two step inexact test
  //reject if bounding box does not intersect bounding box of the disk
  if(!getBoundingBox().intersects(box)) return false;

  //also reject if box does not intersect the plane of the disk
  F32vec4 norm = getNormal();
  F32vec4 half = F32vec4(0.5f);
  F32vec4 boxCenter = half*(box.min+box.max);
  I32vec4 absIntMask = Is32vec4(0x7FFFFFFF,0x7FFFFFFF,0x7FFFFFFF,0x7FFFFFFF);
  F32vec4 absMask = _mm_castsi128_ps(absIntMask);  //used to implement absolute value operation
  F32vec4 absNorm = norm&absMask;
  //maximum change in distance to disk plane over the box (compared to the center)
  F32vec1 maxChange = dotProduct(half*(box.max-box.min),absNorm);
  F32vec1 centerDist = dotProduct(boxCenter,norm) - dotProduct(getCenter(),norm);
  if (_mm_ucomigt_ss((centerDist+maxChange)*(centerDist-maxChange),F32vec1(0))) return false;

  //otherwise we must conservatively assume that they do intersect
  return true;
}

AxisAlignedBox Disk::estimateIntersectionBounds(const AxisAlignedBox &box) const
{
  //return a conservative bounding box by simply intersecting the input bounding box
  //and bounding box of the disk.  This will often not be a tight bounding box for the
  //portion of the disk that lies inside the input box.
  AxisAlignedBox bounds = getBoundingBox();
  bounds.min = simd_max(bounds.min,box.min);
  bounds.max = simd_min(bounds.max,box.max);
  return bounds;
}

void PackedDiskList::setDisk(int index, float centerX, float centerY, float centerZ, 
                                        float normalX, float normalY, float normalZ,
                                        float uAxisX, float uAxisY, float uAxisZ,
                                        float radius, int id)
{
  assert((length > index)&&(index>=0));
  int i = index;
  disks[i].center = F32vec4(0,centerZ,centerY,centerX);
  disks[i].normal = F32vec4(0,normalZ,normalY,normalX);
  disks[i].uAxis = F32vec4(0,uAxisZ,uAxisY,uAxisX);
  disks[i].radius = radius;
  disks[i].id = id;
}

void IntersectPackedDisks(class PrimitiveList *list, class Ray *ray, class IntersectionRecord *hitRec)
{
  F32vec1 tStart = F32vec1(ray->tStart);
  F32vec1 tEnd = F32vec1(ray->tEnd);
  {
    PackedDisk *disks = ((PackedDiskList*)list)->disks;
    for(int count = list->length;count>0;count-=1) {
      F32vec4 norm = disks->normal;
      F32vec4 cent = disks->center;
      F32vec4 uAxis = disks->uAxis;
      //first form the denominator
      F32vec1 denom = rcp_nr(dotProduct(norm,ray->direction));
      //next form the v-axis
      F32vec4 vAxis = crossProduct(norm,uAxis);
      //S is center minus the ray origin
      F32vec4 S = cent - ray->origin;
      //t is the parametric parameter along the ray of the intersection
      //t = (S dot C) / (D dot C)
      F32vec1 t = dotProduct(S,norm)*denom;
      //B is S X D
      F32vec4 B = crossProduct(S,ray->direction);
      //u is the first barycentric coordinate of the intersection
      //u = (vAxis dot B) / (D dot C)
      F32vec1 u = dotProduct(vAxis,B)*denom;
      //v is the second barycentric coorinate of the intersection
      //v = -(uAxis dot B) / (D dot C)
      F32vec1 negV = dotProduct(uAxis,B)*denom;
      //test to see if intersection is valid
      F32vec1 rad2Test = u*u + negV*negV;
      F32vec1 radius2( disks->radius*disks->radius );
      F32vec1 valid = cmple(rad2Test,radius2)&cmpge(t,tStart)&cmple(t,tEnd);
      int test = move_mask(F32vec4(valid))&0x01;
      if (test) { //there is a valid intersection
        F32vec1 invRadius = rcp_nr(F32vec1(disks->radius));
        u *= invRadius;
        negV *= invRadius;
        if (hitRec->func) { //if we have a custom intersection responder, then use it
          (*hitRec->func)(ray,hitRec,F32vec4(t)[0],F32vec4(u)[0],-F32vec4(negV)[0],0,disks->id);
          tEnd = F32vec1(ray->tEnd); //reload in case it has changed
        } else { //otherwise the default is find first intersection
//_asm int 3; //force debugger breakpoint
          hitRec->t = ray->tEnd = F32vec4(t)[0];
          hitRec->u = F32vec4(u)[0];
          hitRec->v = -F32vec4(negV)[0];
          hitRec->w = 0.0f;
          hitRec->primitiveId = disks->id;
          tEnd = t;
        }
      }
      disks++;
    }
  }
}

